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Abstract 

Using Molecular Dynamics simulations based on the effective hamiltonian developed by Zhong, 
Vanderbilt and Rabe [Phys. Rev. Lett. 73, 1861 (1994)] (and fitted on first-principles calculations 
only), the technique of the thermodynamic integration is applied to barium titanate. It allows 
to compute the difference of free energy between macroscopic states with different polarizations, 
from the thermal averages of the forces acting on the local modes. This is achieved by performing 
molecular dynamics under the constraint of fixed polarization. The Landau free energy is thus 
interpreted as a potential of mean force. The thermodynamic integration directly gives access (nu- 
merically) to the Landau free energy of barium titanate as a function of P, without any assumption 
on its analytical form. This technique, mainly used in computational chemistry, allows to make 
a direct connection between phenomenological theories and atomic-scale simulations. It is applied 
and validated to the case of BaTiOa under fixed volume. The results are compared to available 
phenomenological potentials of the litterature. 
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I. INTRODUCTION 



Ferroelectric (FE) solids have been the subject of many theoretical studies for a long time. 
In particular, the nature of their phase transitions, displacive versus order-disorder, and the 
microscopic mechanisms that might be responsible for them, the space and time correlation 
of dipoles, the damping of the soft mode in the high-temperature phases, have attracted 
very soon much interest since these phenomena are of fundamental interest to understand 
the deep nature of ferroelectricity. 

From a phenomenological point of view, FE phase transitions are well described in the 
framework of the Landau theory. This theory assumes the existence of a so-called Landau 
free energy F{N, V, T, P), that can be seen as an "incomplete" thermodynamic potential l||: 

F{N, V, T, P) = -ksTlnZiN, V, T, P), (1) 
in which Z{P) is an incomplete partition function defined by 

Z{N, V, T,P) = J2 e~''/^^^ (2) 
i/p 

In this expression, the summation is extended over the microscopic states (i) with energy 
exhibiting a polarization equal to P. In Landau theory, the equilibrium order parameter P 
is obtained as the value minimizing the Landau free energy. 

The Landau theory has been successfully applied to BaTiOs (BTO) by Devonshire in 
the fiftees . He expanded the Landau free energy in power series up to sixth order 

in the three components of the polarization P^, Py and Pj. Assuming a linear dependence 
with temperature of the coefficient of the P^ term (as stated by Landau theory), he showed 
that the three phase transitions of barium titanate can be described with a relatively simple 
expression of the free energy. Since this pionnering work. Landau theory has been applied 
to barium titanate in many other forms than bulk (filmsjsl, wiresQ]) and to many other 
ferroelectrics 3|- 

The direct calculation of the Landau free energy requires two things: 

(1) performing simulations with fixed polarization, which can be achieved easily within 
constrained Molecular Dynamics (with modified equations of motion obtained from la- 
grangian formalism), or by applying a suitable external electric field 

(2) calculating the thermodynamic potential. 



This last point is the most difficult since standard simulation techniques such as Monte 
Carlo (MC) or Molecular Dynamics (MD) do not give a direct access to the free energy, 
that can not be attained as the thermal average of some microscopic quantity. These two 
techniques only give access to macroscopic states of the system that are minima of the 
thermodynamic potential. 

The coefficients of the free energy expansion for barium titanate have been calculated 

HQ. 

from the first-principles derived effective hamiltonian of Zhong et al[8|, 19| in an elegant way 
by Iniguez et through Monte Carlo simulations under applied external electric field. 
Although Monte Carlo only gives access to the minima of the free energy, these authors used 
the external field to displace the minima in the 3-D space (P^, Py, Pz) and get access to 
the coefficients of the expansion and their temperature dependence. In fact, their approach 
could have been extended to obtain the derivatives of the free energy with respect to the 
polarization at various points, and thus, by integration, get access to the free energy as a 
function of P without any assumption on its analytical form. The approach we propose in 
this article is formally equivalent to this point. 

Indeed, some efficient techniques do exist to get directly insight into the thermodynamic 
potentialicr mo.e precisely into difaencescf thermodynamic potential), such a. umbrella 



sampling 
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free energy perturbation|l2| or thermodynamic integration[l3j. The basic 



idea of thermodynamic integration is to find a "reaction coordinate" A joining continuously 
two macroscopic states Aq and Ai, and calculate the derivative of the free energy with respect 
to A (normally accessible within constrained MD) along this path. This derivative can thus 
be integrated to obtain the difference of free energy between Aq and Ai: 




F(AO-F(Ao)= / i—iX)} d\ (3) 

or even the whole free energy profile along the path. In this formula, F{\) is the incom- 
plete thermodynamic potential defined above. 
In most cases, this yields: 



F(Ai)-F(Ao)= / i^) d\ (4) 

|— I J\a \ / NVT 

in which U is the potential energy [lJ|. 

This technique is widely used in computational chemistry and chemical physics, two 
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fields in which scientists are interested in calculating free energy profiles along a reaction 
coordinate A to obtain the activation free energy of a chemical reaction or its equilibrium 
constant (related to the standart change in free energy) . The free energy for a given value of 
A (referenced to that with Aq) is interpreted as the "potential of mean force" on the reaction 
coordinate 

F(A)-F(Ao) = - f\f{\'))d\', (5) 

i.e. the potential related to the averaged force acting on the reaction coordinate A (usu- 
ally —/(A) is the external force that must be applied to maintain the reaction coordinate 
invariant). 

In our case, if we choose for A the polarization and take as starting point (A=0) the 
paraelectric state of BTO at some given temperature, it is thus formally possible to get 
access to the free energy for any value of P, i.e. the Landau free energy of the ferroelectric 
crystal. Thus we interpret the Landau free energy as a potential of mean force, calculated 
as a function of the polarization, considered as the reaction coordinate. Chemists have been 
using this technique for a long time to calculate reaction free energies and activation free 



energies along paths defined by a reaction coordinate 



3- 



In this paper, we transpose the 



thermodynamic integration technique to the field of ferroelectricity. 



II. THEORETICAL BACKGROUND 



A. Effective Hamiltonian 

We adopt the formalism of the Effective Hamiltonian[^ Sj. In this approximation, the 
number of degrees of freedom of the perovskite (15 per ABO3 unit cell) is reduced to 6 per 
ABO3 unit cell, chosen to be local collective displacements constructed from the lowest- 
energy eigenvectors of the force constant matrix. Two kinds of such modes are considered: 
the "local modes" Ui, that roughly represent the dipole in each unit cell, and the "displace- 
ment modes" Vi, describing the inhomogeneous strain. The use of such degrees of freedom 
has revealed to be efficient and sufficient to describe correctly the thermodynamics of BTO 
(in particular its complex sequence of phase transitions). In terms of these relevant degrees 
of freedom, the potential energy is the Effective Hamiltonian H'^-l'^ {{ui} , {vi} , {ril^}), Tjf 



4 



being the homogeneous strain tensor. 

The local modes Ui are related to the polarization by: 

P=^Y.^^^ (6) 

i 

where f2 is the volume of the supercell and Z* the effective charge of the local modes. 
We define an average local mode u: 

^=^I]^i' (7) 

i 

in which N is the number of unit cells in the supercell. Fixing P is thus equivalent to fix u 
and we use u for convenience in the following (we calculate the free energy as a function as 
u). u has the dimension of a displacement. 



B. Molecular dynamics under fixed polarization 

Our approach in this work is based on molecular dynamics (MD). The first-principles 
derived hamiltonian of Zhong et al can be solved either via Monte Carlo simulations or 
Molecular Dynamics, that should yield similar results under the ergodicity hypothesis. We 
have already validated the MD method with the effective hamiltonian on barium titanatejl^ 
and showed that its complex sequence of phase transitions as well as the temperature evo- 
lution of the polarization and strain can be reproduced successfully as it is in Monte Carlo. 
This is achieved by combining the Nose-Hoover{l3, Q, [3| (to fix the temperature) and 
the Parinello- Rahman [2^ (to fix the pressure/stress) algorithms. More details about this 
molecular dynamics are given in Appendix B. 

Anyway, to calculate differences of free energy along a path indexed by polarization 
(reaction coordinate), it is necessary to perform MD under fixed polarization. Additional 
forces have to be added in the dynamical equations of motion that ensure the constraint 
^ • Ui = Nu. These additional forces can be found in a simple manner starting from the 
lagrangian formalism. Let us call L the lagrangian of the system without constraint, and L' 
that of the system with constraint. 
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in which are Lagrange multipliers and a = x,y and z. Applying the Lagrange equations 
on L' gives a new set of equations of motion for the local modes: 



in which mim is the mass associated to the local modes[16] and the a-component of 
the force acting on the i^^ local mode. Summing over i gives immediately the value of the 
Lagrange multiplier along a: 

i 

Maintaining the polarization fixed is thus equivalent to add to the force acting on each local 
mode (and at each time step) a term equal to (minus) the (space-) average of the forces. The 
total external force that must be applied to maintain P fixed is thus — ff^. 

'Yji fi"^ '^ill be thus qualified in the following of this paper as the (internal) force "acting 
on the polarization" (one must apply its opposite to the system so as to maintain the 
polarization invariant) . 

The equations of motion of the displacement modes Vi are not affected by the constraint. 



C. Thermodynamic integration 

The application of the technique of the thermodynamic integration yields in our case: 

m -Fo = -<fYl (11) 

in which a continuous path is assumed from u = to u. F is the free energy of the entire 
supercell, and the summation is over all the local modes of the supercell. N,V and T are 
assumed constant along this path. A proof of this equation is given in Appendix A. For 
each point of the path, the thermal average of the forces is calculated under fixed u. In the 
previous formula, Fq = F{u = 0). 

In the particular case of a linear path along the a-axis: 

Hua) -Fo = - J2 (^"> (12) 
Jo i 

in which the calculation is performed along the a direction («„ is the a component of u, 
and Ui3^a=0 along the path). 
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The formula can be written in a local form: 

V,F(n) = -5^(/^)(^) (13) 

i 

The Landau free energy is thus the potential of mean force acting on the polarization. It is 
the integral of minus the thermal average of the total force that is acting on the polarization 
{i.e. the total external force that must be added to maintain the constraint). Written as an 
infinitesimal variation 

dF{u) = -Y,{ft){u).du, (14) 

i 

the right hand-side member can be interpreted as a (thermodynamic) work and the link 
with the basic laws of thermodynamics is easy (see Appendix C). 

If we define a free energy per unit cell ^{u), and an average force per unit cell also {Pj, 
we get the very simple expression: 

(/) = -V„^<l> (15) 

We have thus to perform constrained MD simulations at different values of m^, and 
calculate in each case the thermal average of the sum of the forces acting on the local 
modes. This is why MD is well suited to Thermodynamic Integration since the forces are 
directly available. 



III. COMPUTATIONAL DETAILS 



We perform Molecular Dynamics (MD) based on the Effective Hamiltonian of Ref. 
using the code described in Ref. Il6|. We use a 12x12x12 supercell with periodic bound- 
ary conditions for the thermodynamic integration. For the direct MD simulations, we also 
use this size of supercell, excepting below the phase transition, in which finite-size effects 
are responsible for rapid reversal of the polarization: in that case we use 14x14x14 and 
16x16x16 supercells (10 K below the phase transition) to improve the quality of the numer- 
ical results. The equations of motion are solved within the well-known Verlet algorithm. We 



use a Nose-Hoover 



17|,ll8|,ll9fl thermostat to control the temperature. The paths over which 



the integration is performed are sampled with a step of Am = O.OOlao (ao is the BTO lattice 



33l|). Each simulation with fixed polarization is performed 



constant obtained from the LDA 
on 10^ steps, the thermal averages being computed on the 50000 last steps. The time step 



is fixed to 2 X 10~^^ s. All the free energy curves presented herafter are performed at fixed 
strain tensor, and each curve corresponds to a fixed temperature, to ensure we calculate a 
(Helmholtz) free energy curve F{N,V,T,u). 

The thermodynamic integration along the a direction is simply performed, from zero to 
Ua = nAu, via: 

rnAu _ 

/ E > «)^< = ^« E E if'^) (^•^«)' 

i k=0 i 

the thermal average being performed for each value of k under fixed polarization (= kAu). 



IV. THERMODYNAMICS OF BATIO3 UNDER FIXED VOLUME. 

We apply our method to the case of BaTiOs under fixed volume. This system has been 
studied very soon, and it was suggested from Landau theoryj^, that BTO under fixed 
volume might undergo a second-order phase transition at the Curie temperature. Several 



experiments 



21 



22| have confirmed that clamped BaTiOs has the coefficient of (in the 
Landau-Devonshire expansion) positive, making the paraelectric-ferroelectric transition of 
second order. 



A. Isotropic strain 

We first focus on BTO constrained to have a cubic unit cell, with a strain tensor frozen 
to .he va,.es found ,„ ... Mgh-.e,„pe.at„.e cubic phase d.ect MD ..ula.icn.H: 

Vi = V2 = Vs =0.0121 and r]4 = rj^ = ?76=0.0. This corresponds to a lattice constant of 3.995 
A. We note that these strains result from the application of a negative pressure -4.8 GPa (i.e. 
positive stress tensor axx = o'yy = azz =4.8 GPa), which is the one used by Zhong et a/|^ 
to correct the underestimation of the lattice constant by the Local Density Approximation 
(LDA). Moreover, the strain is defined here with respect to this cubic LDA ground state 
(lattice constant 7.46 a.u.). 

Using the Nose-Hoover MD method, we first perform simulations to determine the tem- 
perature evolution of the polarization and the possible phase transitions. In fact we only 
find one phase transition in that case, from cubic (high T) to rhombohedral (low T) (Fig.[T]) 
at ^ 280 K. Zhong et a/fl have already noticed this point through MC simulations on BTO, 
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but with a homogeneous strain tensor frozen entirely to zero. Moreover, the transition seems 
to be second order or at least weakly first-order (the fiuctuations just below the Curie tem- 
perature are very large and finite-size eff'ects make the polarization fiip from a direction to 
another, which makes it difficult to discriminate). 
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FIG. 1: Temperature evolution of the three components of the polarization from direct MD simula- 
tions (without constraint). The polarization minimizing the Landau free energy is shown in orange 
circles (5 calculations surrounding the phase transition). A 12x12x12 supercell is used except 
below the transition where the supercell size is 14x14x14 and 16x16x16 between Tc-IOK and Tc. 

We now turn to thermodynamic integration for a set of selected temperatures surrounding 
the Curie temperature: T=320, 300, 280, 260 and 240 K. We perform the integration along 
two paths starting from u = and going to 0.025aoe2; (along [100]) and to 0.025aQ{e^+ey+ez) 
(along [111]). The free energy, resulting from the integration of the thermal averages of the 
total forces along the second path (i.e. along [111]), is shown on Fig. [2l The minima of the 
free energy are reported on Fig. [T] (orange circles) for the five temperatures: as expected, the 
results obtained by direct Molecular Dynamics correspond to the minima of the free energy. 
We have also computed the profile of the free energy along the [100] direction, leading to 
globally higher free energies than along [111]. 

Interestingly, the rhombohedral ferroelectric phase obtained at low temperature does not 
exhibit shear strain (since shear strain is fixed to zero) as it would be the case for a stress-free 
crystal. Instead, we observe a shear stress, as can be seen on Fig. [5] (inset): the non-diagonal 
components of the stress tensor are not zero at low temperature and decrease to zero at the 
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FIG. 2: Landau free energy of cubic BTO (fixed volume) along [111] for 5 temperatures. FO is 
the free energy of the paraelectric state for each temperature. The five curves plot F{N, V, T, u) — 
F{N, V,T,u = 0) as a function of Ux (= % = Uz) along [111]. 

transition. 



B. Tetragonal strain 

As another example, we freeze the strain to the values found in the tetragonal phase 
at T=280 K, i.e. r]i = 0.01755, r]2 = r]^ = 0.01027 and r]^ = r]^ = r]e = 0.0 (as in the 
previous case, these strains corresponds to a negative hydrostatic pressure of -4.8 GPa at 
this temperature). The corresponding lattice constants are a= 4.017 A and b=c=3.988 A. 
We first perform direct molecular dynamics simulations from T=260 to 400K (Fig. [3|). The 
system slowly evolves from ferroelectric to paraelectric, the transition temperature being at 
^ T=375 K. 

We also calculate under this fixed strain tensor the free energy as a function of (with 
Uy = Uz = 0.0), for eight temperatures (from 260 to 380K) and plot it on Fig. [H The values 
of Ux directly found by molecular dynamics correspond to the minima of the Landau free 
energy curve, as expected (the values minimizing the free energy are reported on Fig. [3l in 
orange circles). Interestingly, the transition is also second order or weakly first-order: within 
the numerical accuracy of our calculations, the evolution of the Landau free energy curves 
as a function of temperature follow the well-known trends of second order transitions, with 
a continuous evolution from a double-well to a single-well (with minimum at zero) curve. 
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FIG. 3: Temperature evolution of the x component of the polarization from direct MD simulations 
(without constraint) in the case of tetragonal BTO (tetragonality along x). The y and z components 
are not shown since they are zero. The polarization minimizing the Landau free energy is shown 
in orange circles (8 calculations). A 12x12x12 supercell is used except below the transition where 
the supercell size is 14x14x14 and 16x16x16 between Tc-IOK and Tc. 
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FIG. 4: Landau free energy of tetragonal BTO (fixed volume) along [100] for 8 temperatures around 
the phase transition. 



In the range of temperature considered (260-400K), we do not find any orthorhombic 



phase. Such a phase natural 
and not under fixed volume 



appears when the simulation is performed under fixed pressure 



16| 



The second-order (or weakly first-order) character of the two phase transitions is illus- 
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trated by the temperature evolution of the stress tensor components (Fig. [5|), that do not 
exhibit discontinuities at the phase transition. 
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FIG. 5: Temperature evolution of the stress tensor in the two cases studied under fixed volume 
(cubic and tetragonal). The strain is frozen in both cases to values obtained from direct MD 
simulations at a negative hydrostatic pressure -4.8 GPa {i.e. positive stress tensor), which explains 
why the stress tensor diagonal components have very high values around 4.8 GPa. 



C. Comparison with phenomenological potentials 

Since the work of Devonshire, the thermodynamics of barium titanate has been the sub- 
ject of many phenomenological studies, and several Landau potentials have been proposed 
as expansions in power series of P up to Q^^ or 8*'' order for this material. The Landau 
potential is more often given as a Gibbs free energy and for a stress-free crystal {cr=0). 



Converting to an Helmholtz free energy 
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2J|, F - Fq{= AF) usually writes as the 



sum of three terms, in which it appears as a function of temperature, strain tensor ry and 
polarization: 

AF{rj, T- P) = Fl{T; P) + F,i{r^, T) + F,(r/, T; P), (16) 
in which rj denotes the strain tensor (given in Voigt notations) and 
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F,(T; P) = a,{Pl + + P|) + a^P'^ + P^ + P,^) 
+aUP'Py + P'P' + PyP!) + o^iiiiP! + i^,' + n') 
+o.ii2{Pl{Pt + + ^.'(^' + Pt) + + Pt)) 

+au3P!p;p! + + PS + n') 

+ani2(P.'(P,' + P|) + P,'(PJ + P|) + PUP^ + P;)) 

+0^ll22{PX + P'Pt + P'yPt) 
+an23iP^PSP! + Py'P^P! + PtPlP^) 



+Ci2(r/i?72 + ViVs + V2V3) + ^ivl + vi + Ve), 



where the Cij are the elastic constants, 



F,{ri, T- P) = -qii{r^iPl + r^a^,' + ^Pl) 

-quimiP; + P^) + miPl + P'.) + n^iPl + P'y)) 

-q^iVePxPy + V^PxPz + TliPyPz), 



where the qij are the electrostrictive coefficients. 

The coefficients entering the expression of are the same as those entering the 
free energy of the stress-free crystal, excepting the 4*'^-order coefficients, that must be 
renormalized{23, [2^. These renormalized 4*''-order coefficients al-^ and are calculated 
from ail and ai2 (from the stress-free crystal) and from the elastic and electrostrictive 



23l to calculate these renor- 



coefficients[23||. In the following, we will use the values of Ref. 
malized coefficients. 

This function F — Fq should be directly comparable to the free energy we have calculated 
by thermodynamic integration. We will restrain to the comparison with two recent studies 
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that have proposed 8*'*-order expansions of the Landau free energy: Li et a/[25|| and Wang 



et al 
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27|. However, the reader must be cautionned that the comparison is not always 
relevant for the following reasons: 

(i) the strains have not the same reference in the two models: in the effective hamiltonian 
the strain is defined with respect to the cubic Pm3m LDA ground state (lattice constant 
7.46 a.u.), which causes problems in the comparison of the quadratic coefficients of F. 

(ii) the effective hamiltonian has a Curie temperature ~ 120 K lower than the true one. 
Thus the quadratic coefficient, that is expected to evolve with T as a(T — To), will be 
naturally strongly shifted with respect to the phenomenological potentials, that reproduce 
much better the transition temperatures. 

(iii) finally, some free energy curves can give very different coefficients according to the 
degree of the polynomial function with which they are fitted. 

For these reasons, we have prefered to plot the quadratic coefficient as obtained from our 
calculations without correction. From a general point of view, the coefficients fitted on our 
data have orders of magnitude in good agreement with the phenomenological potentials, 
although their temperature evolution appears as more complex, as already pointed out by 
Iniguez et al[l|. 

1. Isotropic strain case 

In the first case studied (isotropic strain), we have = Py = Pz and ?7i = = Thus 
the free energy writes: 



AF(r7, T- P,) = 3(ai - (gn + '2qu)Vi)P! 
+3(a?i + a?2)^x + 3(aiii + 2au2 + ai2z/^)P^ 
+3(aiiii + 2aiii2 + 01122 + "1123)^^^ + 



The curves of Fig. [2] have been fitted on S*'' order polynomial functions (without odd 
term) . 

The coefficients of the high-order terms give access to complex combinations of the a 
coefficients (aim + '^ctiiu + ^1122 + ttii23 and am + 1a\x2 + cti23/3). The coefficient of P^ 
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is plotted on Fig. [6]-a together with + calculated from the potentials of Wang et al 
and Li et al. It falls just between the two and has the same sign as that of Wang. 
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FIG. 6: Temperature evolution of al^ + (^) ^i^d a^^ (b) according to Wang et al (orange 
diamonds) and Li et al (blue circles). Red diamonds: coefficient of P^, calculated from the isotropic 
strain case and divided by 3 (a), and calculated from the tetragonal strain case (b). 

The coefficient of is plotted on Fig. [Tl-a as well as ai calculated from the potentials of 
Wang et al and Li et al. As already said, they are not directly comparable since under fixed 
strain, ai should be renormalized by a constant term ai —rji^qn + 2^12). Anyway the slopes 
are close. The term found from our calculation approximately varies with T as a{T — Tq) 
with To ~ 270K, which roughly corresponds to the transition temperature in the isotropic 
strain case. 



2. Tetragonal case 

In the tetragonal case, we have Py = Pz = and ri2 = rj^. Thus, the Landau free energy 
writes: 
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FIG. 7: Temperature evolution of a\ (a and b) according to Wang et al (orange diamonds) and Li 
et al (blue circles). Red diamonds: coefficient of P^, calculated from the isotropic strain case (a), 
and calculated from the tetragonal strain case (b). 



AF{r], T- P,.) = aiP^ + a^.P^ + aniP^ + aniiP^ 
-qnViP! - '2qi2mP^ + Fei 



Then, as in the previous section, v^^e have fitted the curves of Fig. [4] obtained by the 
thermodynamic integration by polynomial functions up to 8*^ order. For each temperature, 
we get a set of 4 coefficients that we now compare to the ones of Wang el al and Li et al. 

The highest order terms should directly correspond to the am and aim coefficients. The 
values provided by our data correspond very well to the phenomenological ones above ~ 325 
K (Fig. [8]). Below this temperature, these coefficients strongly vary. 

The variation of the P^ coefficient, that should compare directly to a^^, reveals more 
complex (Fig. [H-b), as well as that of the quadratic one (Fig. [Tl-b). 
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FIG. 8: Temperature evolution of the ami and am coefficients according to Wang et al (orange 
diamonds) and Li et al (blue circles). Red diamonds: coefficients of the (a) and (b) terms 
as obtained from our calculations in the tetragonal strain case. 

V. CONCLUSION 



In this work, we have proposed and validated a method to calculate numerically the 
Landau free energy of a ferroelectric material, without any assumption on its analytical 
form and using the effective hamiltonian of Zhong et a/j^ directly fitted on first-principles 
calculations. It is based on the technique of the thermodynamic integration coupled with 
Molecular Dynamics performed at fixed temperature (Nose-Hoover) and fixed polarization 
along paths in the 3-D space (P^, Pj^, P^). Although possible also with Monte Carlo sim- 
ulations. Molecular Dynamics is well suited to this kind of problem since thermodynamic 
integration requires to calculate the forces, directly available within MD. 

The method has been applied to the case of BTO under fixed volume and showed that 
minima of the calculated Landau free energy correspond to the results of direct MD simu- 
lations. The results have been compared with some available phenomenological potentials. 
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This method connects directly the phenomenological theories to the atomic-scale simu- 
lations. The next step is to extend the simulations to the isothermal-isobaric ensemble to 
obtain the Gibbs free energy as a function of polarization. This will be achieved by coupling 
Nose-Hoover, Parinello-Rahman and also fixed polarization molecular dynamics. Other ex- 
tensions of this approach can be imagined, such as the computation of entropy profiles as a 
function of P, the application to low-dimensional systems and to complex perovskite (with 
antiferrodistortions) and multiferroic materials, where the free energy could be computed as 
a function of polarization and tilt angle or magnetization j28|. 

The method also offers a wide freedom to use any internal parameter of the system as 
variable to expand the free energy 



VI. APPENDIX A: THE LANDAU FREE ENERGY AS A POTENTIAL OF MEAN 
FORCE 

In this section we give proofs for Eq. [HI very similar to demonstrations that can be found 

nnnn 

in the chemical scientific litterature on thermodynamic integration[13|, HJ, 115|, 129||. Let us 
considere two macroscopic states of the system (1) and (2), each defined by (N,V,T), but 
having different polarizations P = Z* /QJ^i'^i = NZ*/Qu. It is more convenient to work 
with the average local mode u, related to P only by a constant factor. 

We introduce as usual an incomplete free energy, for a fixed value of u, defined by 

F{N, V, T, u) = -kBTlnZ{N, V, T, u), (17) 

Z{N,V,T,u) being an incomplete partition function defined by a summation over the mi- 
croscopic states having a polarization equal to P = NZ* /Qu. 
The difference of free energy between state (1) and state (2) is 

F(N, V, T, U2) - F{N, V,T,u,)= f .du, (18) 

—t 

in which the derivative with respect to a vector denotes the gradient operator = 
{d/dux, d/duy, d/ duz) . § is the integral over a continuous path joining (1) to (2) in the 
3-D space {u^, Uy, u^). 

We have thus to calculate the term 

dF ksT dZ 



du Z du 
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(19) 



We start by calculating the incomplete partition function Z{N, V,T,u) 

Z{N, V, T,u) = J2 e"^'^'' (20) 

r/u 

(r) denoting the microscopic states and the summation being over the (r) with polar- 
ization P = NZ*/Qu and /? = l/ksT. We use the classical approximation to calculate Z 
and omit in the following N,V,T in the notations, all the partial derivatives being performed 
at constant V, T. In this framework, the partition function is an integral over the phase 
space, the summation being restricted to the microscopic states satisfying the constraint 
^^Ui = Nu. It writes fllaO]: 



^(^) ~ TW / dpi---dpN I dui...duM 



I 



In this expression, the energy of the microstates is: 

N 



E{pi...pN, Ui...Un) = + H^^^{Ui...Un) 



i=l 

The Pi are the conjuguate momenta of the local modes Wj. In the expression of N\ 
is not introduced since the local modes, which correspond to well identified unit cells, are 
discernable. Moreover we do not write the displacement modes Vi for simplicity. 

Integrating over the pi yields: 

f 2TT'mkRT 1 '^^^^ 



^(^) = |~~~^^^| J dui...du 
X(5(^ Ui — Nu)e 



N 

l3H''ff{ui...UN) 



This integral can be calculated in several manners. One can remove the Dirac function 
by replacing one of the Ui by Nu — YljjLi^j- For instance. 



f 2'KmkBT 



Z{u) = < — > / dui...du 



3N/2 



I'M-! 



xe' 
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The microscopic states satisfying the constraint thus appear as functions of u and in- 
dexed by ui,U2, ■■■Un-i,Pi---Pn- This is more obvious by rewriting: Ht^J Ju) — 



f 2nmkBT 



'^('") ^ I I / dui...du 

Xe «l--"Ar-i^ ' 



3N/2 



The derivative with respect to u is: 

OF ksTdZ 



du Z du 

dHl^^ (u) 

dui...duN-i 

ou 

-PHZff ^ (u) 

/ dui...duN-ie «i-"jv-i^ ' 



X 



which can be rewritten: 



'd^^JfiN / dpi...dpN / dui...duNS(^Ui- Nu) 

i 

X TT^ ^ 

du Z 

One recognizes the thermal average (under fixed polarization) in the canonical ensemble 

dHtf ^ - (u) 

of the microscopic quantity — "^"^5"^ — . 
This quantity is in fact: 

ou OUn 

in which /n is the force acting on the N*'* local mode. 
And finally, 

^--N{h) (22) 

It is clear that the N*'* local mode does not have a particular role. It appears here because 
of the choice made to calculate the initial integral. If we had chosen to remove Ui instead 
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of Mat, we would have found 

In fact, the system is completely homogeneous and the thermal average of the forces does 
not depend on the site on which it is calculated: 

/i) = (/I) = - = (/^) (24) 

To improve the statistics of the thermal average, it is more interesting to calculate the 
thermal averages of all the forces and to sum them: 

dF\ 



du I 



) N,V,T «=1 



The thermal average in this expression is performed at fixed polarization. The quantity 
that must be integrated is thus the mean force acting on the "reaction coordinate", i.e. 
the polarization (this is the total force that has to be added in the equations of motion to 
maintain the polarization constant) . The free energy is thus the potential of the mean force 
acting on the polarization. 

We note that in some cases, the use of more complex reaction coordinates can give rise 
to more complex expressions for the derivative of the free energy jl^. [s^. This is not the 
case here. 

Very similar arguments in the isothermal-isobaric ensemble lead to 

N 

E {f^ ' (26) 

N,P,T «=1 

which is formally the same expression, except that the thermal average must be performed 
under fixed temperature, pressure and u. This simple extension is related to the fact that 
local modes and homogeneous strain are independent variables of the Effective Hamiltonian. 

In term of polarization (instead of average local mode), it writes: 





NZ 

N,P,T «=1 



The calculation of |^ could be performed by using the Parinello-Rahmanj2^ scheme 
combined to Nose-Hoover and fixed polarization. 
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VII. APPENDIX B: MOLECULAR DYNAMICS USING EFFECTIVE HAMILTO- 
NIANS 

The effective hamiltonian is solved in our approacfi by Molecular Dynamics [l^. We 
describe here the main features of this approach. A specific mass is affected to the local 
modes and displacement modes, and these degrees of freedom evolve in time according to 
the Newton equations of motion. At each time, the forces acting on these local modes derive 
from the effective hamiltonian: 

where Im and dsp stand respectively for local mode and displacement mode. These 



equations are integrated numerically within the well-known Verlet algorithm [3jJ: the local 
mode at time t + h (h is the time step) is computed from the local modes at t and t — h, 
using the forces calculated at t: 

u,{t + h) = 2u,{t) -u,{t~h) + ^fl'^it) (30) 
and idem for the displacement modes. 

To fix the temperature, we have implemented the Nose-Hoover method 



17|,Q,Q|. This 



method consists in adding to the system a fictitious degree of freedom s representing the 
thermostat, and that couples to the equations of motion, according to: 

d'^Vi TtJon i / .N dvi 



ft^ -rridspCit)-^ 



dt^ "''^''dt 



dt Q ^T{t) 



where T(t) and Tq are respectively the instantaneous and targetted temperature, g is the 

dint 
dt 



number of degrees of freedom in the system, and C,{t) = ^^(t). Q is the so-called Nose 
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mass, that must be chosen correctly so that the exchange of energy between the system and 
the thermostat are neither too slow nor too fast. The velocities at time step t are calculated 
from the Ferrario algorithm {32 1: 

duj _ 3ui(t) - Aujjt - h) +Ui(t- 2h) 

~dt^^^ ~ 2h ^"^^^ 

VIII. APPENDIX C: RELATIONSHIP TO THERMODYNAMICS 

The basic laws of thermodynamics allow to write: 

dF = 6W' - PdV - SdT (32) 

where 5W' is the work of external forces other than pressure forces, and —PdV is the 
work of pressure forces. A system evolving reversibly under constant N,V and T obeys: 
dF = 6W'. 

In our case, the variation of the mean local mode from m to u + du (which has the 
dimension of a distance) is related to the work of the force — J2i {ji^^-, ^-S-- 

dF = 5W' = -Y,{fl''')-du (33) 

i 

The Landau free energy variation is thus the work of (minus) the mean force acting on 
the polarization (this is the external force added to maintain the constraint in the molecular 
dynamics). 

For the Gibbs free energy: 

dG = SW + VdP - SdT (34) 

and a system evolving reversibly under constant N,P and T obeys: dG = 6W'. Thus in such 
conditions 

dG = 6W' = -J2 (fl"") -du (35) 



[1] J. Iniguez, S. Ivantchev, J. M. Perez-Mato, A. Garcia, Phys. Rev. B 63 (2001), 144103. 



23 



[2] A. F. Devonshire, Phil. Mag. 40, 1040 (1949). 
[3] A. F. Devonshire, Phil. Mag. 42, 1065 (1951). 
[4] A. F. Devonshire, Adv. Phys. 3, 85 (1954). 

[5] V. B. Shirokov, Y. I. Yuzyuk, B. Dkhil, V. V. Lemanov, Phys. Rev. B 75, 224116 (2007). 
[6] J. Hong, D. Fang, Appl. Phys. Lett. 92, 012906 (2008). 

[7] "Physics of Ferroelectrics, a modern perspective", K. M. Rabe, C. H. Ahn, J.-M. Triscone Eds, 
in Topics in Applied Physics, volume 105. Chapter "A Landau primer for ferroelectrics", P. 
Chandra, P. B. Littlewood. See also the Appendix A by L. Q. Chen. 
[8] W. Zheng, D. Vanderbilt, K. Rabe, Phys. Rev. Lett. 73, 1861 (1994). 
[9] W. Zhong, D. Vanderbilt, K. Rabe, Phys. Rev. B 52 (1995), 6301. 
[10] L. Li, L Vorobyov, T. W. Allen, J. Phys. Chem. B 112 (2008),9574. 
[11] G. M. Torrie, J. P. Valleau, J. Comput. Phys. 23 (1977), 187. 
[12] R. W. Zwanzig, J. Chem. Phys. 22 (1954), 1420-1426. 

[13] J. G. Kirkwood, in "Theory of Liquids", edited by B. J. Alder (Gordon and Breach, New York, 
1968). 

[14] E. Darve, A. Pohorille, J. Chem. Phys. 115 (2001), 9169. 

[15] See for instance: T. P. Straatsma, H. J. C. Berendsen, J. Chem. Phys. 89 (1988), 5876; N. F. 

A. Van Der \^gt, W. J. Briels, J. Chem. Phys. 109 (1998), 7578; P. A. Apte, L Kusaka, Phys. 

Rev. E 73, 016704 (2006). 
[16] G. Geneste, Comput. Mater. Sci., submitted. 
[17] S. Nose, J. Chem. Phys. 81 (1984), 511. 
[18] S. Nose, Mol. Phys. 57 (1986), 187. 
[19] W. Hoover, Phys. Rev. A 31 (1985), 1695. 
[20] M. Parinello, A. Rahman, Phys. Rev. Lett. 45 (1980), 1196. 
[21] M. E. Drougard, R. Landauer, D. R. Young, Phys. Rev. 98 (1955), 1010. 
[22] E. Stern, A. Lurio, Phys. Rev. 123 (1961), 117. 
[23] J. Hlinka, P. Marton, Phys. Rev. B 74 (2006), 104104. 
[24] S. Nambu, D. A. Sagala, Phys. Rev. B 50 (1994), 5838. 
[25] Y. L. Li, L. E. Cross, L. Q. Chen, J. Appl. Phys. 98 (2005), 064101. 

[26] Y. L. Wang, A. K. Tagantsev, D. Damjanovic, N. Setter, V. K. Yarmarkin, A. L Sokolov, 
Phys. Rev. B 73 (2005), 132103. 

24 



[27] Y. L. Wang, A. K. Tagantsev, D. Damjanovic, N. Setter, V. K. Yarmarkin, A. I. Sokolov, J. 

Appl. Phys. 101 (2007), 104115. 
[28] I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, L. Bellaiche, Phys. Rev. Lett. 99 (2007), 

227602. 

[29] M. Mezei, J. Chem. Phys. 86 (1987), 7084. 

[30] W. K. Den Otter, J. Chem. Phys. 112 (2000), 7283. 

[31] L. Verlet, Phys. Rev. 159 (1967), 98. 

[32] M. Ferrario, J. P. Ryckaert, Mol. Phys. 54 (1985), 587. 

[33] ao = 7.46 a.u. from Ref. [9 



25 



